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! Abstract 

00 

We discuss the choice of parameters and report some results for 
unquenched simulations of the Schrodinger functional with a non- 
^3 '. hermitean variant of Liischer's multi-boson algorithm. 

Today Hybrid Monte Carlo (HMC) is the standard algorithm employed 
! for simulations with dynamical fermions. In spite of its general success it 

seems desirable to have other methods at one's disposal. In particular the 
multi-boson technique proposed by Liischer seems interesting. Apart from 
• its theoretical appeal one may perform consistency checks and hope for bet- 

ter efficiency, in particular with regard to slow topological modes [01. Better 
^ , numerical stability and more flexibility in the treatment of statistical prob- 

^ \ lems with exceptional configurations may be further advantages. Soon 

after Liischer's proposal a non-hermitean variant of the algorithm has been 
advocated and initial tests have been performed 0, which we extend here. 
Experiments with the original proposal are reported in Pj. 

The contribution to the QCD Boltzmann factor from two flavors of dy- 
namical quarks is given by 

exp[-5q,arks(f/)] = det(M)2 = I det{M)\\ (1) 

where M = M{U) is the (sofar unimproved) Wilson Dirac operator with 
the hermiticity property M"!" = 75M75. For the multi-boson algorithm we 
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employ a polynomial P{M) which, over the spectrum of M, approximates 
the inverse, 

|det(M)p = |detP(M)|-Vet(l-^(M))|^ (2) 

such that R is a small remainder. This enables us to represent the dominant 
part of ^quarks a bosonic path integral 

|detP(M)|-2 _ 

J D0D0texp|-5]|(M-Zfc)0fe(x)|4, (3) 

I kx J 

where Zk, k = 1, . . . ,n, are the roots of P. We now update by a sequence of 

• some proposal {U, 0) {If, 0') obeying detailed balance with respect 
to the sum of (Q) and the gluon action 

• acceptance with probability q{R, R', x), 

where R, R' are the old and new remainders and x is a complex random field 
governed by some probability distribution p{x)^ i^i our case p oc exp(— x'I'x). 
This compound can be proved to be a valid algorithm if 

j#§4^ = |det(»')|-^ (4) 

{q{R\R, x))x 

with W = {1 — R')^^{1 — R) holds stochastically (i.e. averaged over A 
simple (non-stochastic) solution would be 

go = min(l, |detiy|-2). (5) 

It requires the computation of the det of W. We here use the "noisy algo- 
rithm" of corresponding to 

To evaluate q the apphcation of W to vectors suffices, and the required 
inversion of 1 — i?' with some inverter like BiCGstab is rather uncritical. For 
completeness we mention that the variant called "non-noisy" in was found 
incorrect in the implementation described there. 
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Following we construct P by using Chebysliev polynomials for R. On 
families of nested ellipses with centers at d, 

z {9,ip) = d — ecosh.{9 + iip), (7) 

they approximate the inverse with a rate \R\ < exp(— + l)(6'o — 9)). 
Here d, e are fixed parameters, 9 labels the ellipses and (p G [0, 27r) traces 
them. The polynomial is determined (up to a factor) by the roots Zk = 
z{9o, 2'Kk/{n + 1)), A; = 1,2, ... ,n, lying on the ellipse passing through the 
origin, ecosh^o = d- Due to even-odd symmetry, the spectrum of M is 
symmetric under A ^ 2 — A and we hence set d = 1. 

To implement the correction step we have to evaluate Rx- The factorized 
form, 1 — R = CnMY[{M — Zk), tends to be numerically unstable 0, 0. 
Here it can be avoided and replaced by a uniformly stable two step recursion 
starting from Xo = X leading to Xn+i = Rx- ^ is straightforwardly based 
on the standard recurrence for Chebyshev polynomials. The intermediate Xk 
are the remainders for lower degree polynomials. We follow the recursion 
to investigate the choice of degree n and the focal distance e. With some 
trial parameters we produced some equilibrated [/-configurations, and for 
one of them Figs.|I],^ show the quality of approximation. We see that 
asymptotically the best inversion is achieved for e ~ 0.6 which implies an 
oblate spectrum. For Monte Carlo application, however, n ~ 20 turned 
out to lead to about optimal results. In this range the value of e is rather 
uncritical, and e = 0, where the ellipses degenerate to circles (ecosh6' held 
fixed) and the polynomial to the geometric series, is an acceptable choice. 
This is also confirmed by some direct Monte Carlo runs. In summary, we 
found it practical to use inversion as a tool to infer the spectral information 
necessary to determine the parameters for simulation. The emerging picture 
was stable for various gauge fields and random x that were tested. 

Under even-odd preconditioning one replaces M by M = 1 — MoeM^o 
with det M = det M, and M^e, M^o are blocks of M connecting the even/odd 
sublattices. An application of M has the same complexity as M, but it is 
better conditioned. A pair of eigenvalues A, 2 — A of M is mapped on one 
eigenvalue of M given by their product. Under this mapping ellipses with 
parameters d = l,e are mapped to ellipses with d = 1 — e^/2, e = e^/2. In 
this way the optimal parameters for inversion of M are given by the e optimal 
for M. It again turns out, for the lattice parameters of Figs.|l|,0, that for the 
n relevant for efficient simulation, e = is close to optimal. The errors for 
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Figure 1: Remainder of the inversion with e = for the Schrodinger func- 
tional, L = T = 4, /5 = 6.4, K = 0.15,9 = tt/S, background field "A" [9] 

both cases (same degree n) are connected as |-R| ~ l-Rpj which implies a 
much improved approximation for M. It is interesting that one can prove 
the relations (for even n) 

det{M - Zk) = det{M - Zk') (8) 
det(l - ^) = const X det(l - i?) (9) 

where k' is some permutation of k. Although R is much smaller than R, 
we get (for every single [/-field) the same weight from the boson fields. As 
observed in this allows us to stick to M for the boson fields, which yields 
a much simpler structure for their influence on [/-updating. Would we use 
(H) for the acceptance step, then also this would be identical for R or R. 
In the stochastic case with (^, however, preconditioning dramatically raises 
the acceptance such that n may effectively be about halved. This is due 
to reduced fluctuations in R as compared to R. It is trivial to derive the 
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Figure 2: As Fig.l, ratio of remainders for several e- values 
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{q{R,R',x))x<{qo{R,R')) 



(10) 



and its preconditioned analog. One may thus estimate the loss in acceptance 
from the noisy method which turned out to be tolerable in the preconditioned 
data given below (49% down from 75% with a.t n = 8). 

In Fig.^ results of several multi-boson simulations are shown together 
with a result from preconditioned HMC for the same parameter s^]. They 
are obviously completely consistent for a range of acceptances with and with- 
out preconditioning. The autocorrelation times are given in the table. All 
r refer to units of 1000 M0 applications. The proposals are generated with 
a certain combination of microcanonical and heatbath sweeps. While the 
multi-boson algorithm seems advantageous for the plaquette, there is an ad- 
vantage to HMC for the eigenvalue. In actual CPU time on the Quadrics Ql 
the multi-boson algorithm is faster for both quantities for our particular im- 



5 



plaqucttc 


n 


eigenvalue 




24 






(no pre-cdx) 






6 


X. 


X 


8 

12 


X 


u 

III 


HMC 


L 

1 ' ' 1 



0.661 0.662 0.078 0.08 



Figure 3: Results for plaquette and lowest eigenvalue of M'^M 
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plementations. We plan to clarify this issue further by a simulation on an 8^ 
lattice, but it seems likely that without further new ideas there are no large 
factors in efficiency attainable between the two rather different methods. 
I would like to thank Burkhard Bunk for discussions. 
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